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ABSTRACT 



Context. The radiative energy balance in the solar chromosphere is dominated by strong spectral lines that are formed out of LTE. It 
is computationally prohibitive to solve the full equations of radiative transfer and statistical equilibrium in 3D time dependent MHD 
simulations. 

Aims. To find simple recipes to compute the radiative energy balance in the dominant lines under solar chromospheric conditions. 
Methods. We use detailed calculations in time-dependent and 2D MHD snapshots to derive empirical formulae for the radiative 
cooling and heating. 

Results. The radiative cooling in neutral hydrogen lines and the Lyman continuum, the H and K and intrared triplet lines of singly 
ionized calcium and the h and k lines of singly ionized magnesium can be written as a product of an optically thin emission (dependent 
on temperature), an escape probability (dependent on column mass) and an ionization fraction (dependent on temperature). In the cool 
pockets of the chromosphere the same transitions contribute to the heating of the gas and similar formulae can be derived for these 
processes. We finally derive a simple recipe for the radiative heating of the chromosphere from incoming coronal radiation. We 
compare our recipes with the detailed results and comment on the accuracy and applicability of the recipes. 

Key words, radiative transfer - Sun: chromosphere 
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1. Introduction 

The chromospheric radiative energy balance is dominated by 
a small number of strong lines from neutral h ydrogen, singly 
ionized calcium and singly ionized magnesium (IVernazza et alJ 
although a large number of iron li nes may be equall y 
important in the lower atmosphere ( A nderson & Ath ay 1989). 
The strong lines are formed out of Local Thermodynamic 
Equilibrium (LTE) which means that numerical simulations 
of the dynamic chromosphere need to solve simultane- 
ously the (magneto-) hydrodynamic equations, radiative trans- 
fer equations and rate equations for all radiative transi- 
tions and energy levels involved. This is a major compu- 
tational eflTort that has been accomplished in one dimen- 
sional simulatio ns of w ave propagation in the quiet solar 
atmosphere (Klei n et al l' 1976, 1978; Carlsson & Stein 1992, 
1994, 1995, 1997, 2002; Ramma cher & Ulmsc hneider 2003; 
Fossum & Carls son 2005, 2006), sunspots dBard & Carlsson. 
20101). solar flar es ( Abbett & Hawley 1999: lAllred et al.ll2005l: 



Kasparova et al. j2009l: iCheng et al.i i2010i) and stellar flares 



Kasparova et ai. jzi 
dAllred et al. 200^ 



In three-dimensional simulations the computational work re- 
quired to solve the full set of equations becomes prohibitive; 
there is a need for a simplified description that retains as much as 
possible of the basic physics but brings the computational work 
to tractable levels. The aim of this paper is to provide such sim- 
plified recipes for the chromospheric radiative energy balance. 

In addition to the cooling in strong lines and continua the 
same transitions may provide heating of cool pockets in the chro- 
mosphere. We also get heating in the chromosphere from coro- 
nal radiation impinging from above. We will here address all 



three processes: cooling from lines and continua, heating from 
the same transitions and heating through absorption of coronal 
radiation. 

We use detailed radiative transfer calculations from dynami- 
cal snapshots to derive simple lookup tables that only use easily 
obtainable quantities from a single snapshot from a simulation. 
We test these recipes by applying them to several simulation runs 
and compare the radiative losses from detailed solutions with 
those from the recipes and find in general good agreement. 

The structure of the paper is as follows: In Sec. [2] we describe 
the hydrodynamical snapshots used for deriving the lookup ta- 
bles, in Sec. [3] we describe the atomic models used, in Sec.lHwe 
derive the recipes for radiative cooling from strong lines of hy- 
drogen, calcium and magnesium, in Sec. [5] we derive the recipes 
for heating of cool pockets of gas from the same lines, in Sec.[6l 
we compare the recipes with a detailed calculation, in Sec. [7] we 
derive the recipes for the heating by incident coronal radiation 
and in Sec. [8] we end with our conclusions. 



2. Atmospheric models 

To calculate the semi-empirical recipes we use two sets of hy- 
drodynamical snapshots. The first is a simulation of wave propa- 
gation in the solar atmosphere using the RADYN code (Carlsson 
& Stein 1992, 1994, 1995, 1997, 2002). RADYN solves the one- 
dimensional equations of conservation of mass, momentum, en- 
ergy and charge together with the non-LTE radiative transfer 
and population rate equations, implicitly on an adaptive mesh. 
The grid positions of the adaptive mesh are determined by a grid 
equation that specifies where high resolution is needed (typically 
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where there are large gradients). The grid equation is solved si- 
multaneously with the other equations. The simulation includes 
a 5 -level plus continuum hydrogen atom, 5 -levels of singly ion- 
ized calcium and one level of doubly ionized calcium and a 
9-level helium atom including 5 levels of neutral helium, 3 of 
singly ionized helium and one level of doubly ionized helium. 
All allowed radiative transitions between these levels are in- 
cluded. At the bottom boundary we introduce a velocity field that 
was determined from MDI observations of the Doppler shifts of 
the 676.78 nm Ni i line. Coronal irradiati on of the chromosphere 
was included us ing the prescri ption of IWahlstrom & CarlssonI 
(1994) based on Tobiskal (11991 ). The eff'ects of non-equilibrium 
ionization are handled se lf-consistent l y. The simulation is es- 
sentially the one used in lJudge et aP (120031) but with updated 
atomic data for calcium. A number of ingredients that are poten- 
tially important for the chromospheric energy balance are miss- 
ing: line blanketing is not included (especially the iron lines may 
be important ( An derson & Athayl [l 9 89) ) . cooling in the strong 
Mg II h and k lines is not included and the description of partial 
redistribution is very crude for hydrogen and altogether missing 
for calcium. Most important is probably the restriction of only 
one dimension. This prevents shocks from spreading out, forces 
shock-merging of overtaking shocks and prevents the very low 
temperatures found in 2D and 3D simulations from strong ex- 
pansion in more than one dimension. 



The other snapshot is from a 2D simulation of the solar atmo- 
sphere spanning from the upper convection zone to the corona. 
This simulation was perfor med with the radiation-magneto- 
hydrodynamics code Bi frost (iGudiksen et al . 201 1) and includes 
the eff'ects of non-equilibrium ionization of hydrogen using the 
methods described in Leenaarts et al.. (2007). A detailed de- 
sc ription of the s e tup an d results of this simulation is given 
in iLeenaarts et al ] (l2011h . Also here we have shortcomings in 
the description. Line blanketing is included but with a two- 
level formulation for the scattering probability, hydrogen non- 
equilibrium ionization is included but the excitation balance 
(needed to calculate the individual lines) is very crude. 



Both simulations thus have their shortcomings. The purpose 
of this paper is, however, to develop recipes that reproduce the 
results of detailed non-LTE calculations. It is thus not crucial that 
the simulations reproduce solar observations in detail as long as 
they cover the parameter space expected in the solar chromo- 
sphere. 

For hydrogen it is important to include the coupling be- 
tween the non-equilibrium ionization and the non-LTE excita- 
tion and we therefore use the RADYN simulation for the hydro- 
gen recipes. For magnesium and calcium it is more important 
to include PRD and have a large temperature range than to in- 
clude eff'ects of non-equilibrium ionization and we choose the 
Bi frost 2D simulatio ns. Non-equilibrium ionization is n ot im- 
portant for calcium ( Wedemever-Bohm & C arlssonll201 ll) while 
iRammacher & Ulmschneider ( 2003) and test-calculations with 
RADYN indicate that ionization/recombination times are as long 
as 100-500 s for magnesium and therefore of importance for the 
ionization balance. However, we will see that magnesium is al- 
most exclusively in the form of Mg ii in the temperature range 
where magnesium is important for the cooling/heating of the 
chromosphere and the neglect of non-equilibrium ionization is 
therefore not important for the results. 



CA II 




Fig. 1. Term diagram of the Ca ii model atom. The continuum is 
not shown. Energy levels (horizontal lines), bound-bound tran- 
sitions (solid) and bound-free transitions (dashed). 



3. Atomic models 

3.1. Hydrogen 

We take the population densities directly from the RADYN 
simulation and thus use the same atomic parameters. We use 
a five-level plus continuum model atom with energies from 
Bashkin & St oned (Il975h . oscillator strengths from JohnsoiJ 
(1972), photoionization cross-sections from Menzel & PekerisI 
(193 51) and coUi s ional excitation and ionization rates from 
Vriens & Smeets (1980) except for Lyman-a where we use 
Janev et al. (1987). We use Voigt profiles with complete redis- 
tribution (CRD) for all lines. For the Lyman lines, profiles trun- 
cated at six Doppler widths are used to mimic eff'ects of partial 
redistribution (PRD) (iMilkev & Mihalaslll973h . 



3.2. Calcium 

We use a five-level plus continuum model atom. It contains 
the five lowest energy levels of Can and the Cam contin- 
uum (see Fig. [T]). The energies are fro n i the NIST database 
with data coming from 'Sugar & CorlissI (Il985l) . the oscillator 
strengths are from Theodosiou (1989) with lifet i mes a gree- 
ing well with the experiments of Jin & ChurchI (Il993h and 
the broadening par ameters are from the Vienna Atoniic Lin e 
Database (VALD) (iPiskunov et al.1 119951: iKupka et al.1 Il999l) . 
Photoionization cross-sections are fro m the Opa city project. 
CoUisional excitation rates are from Mel endez et al. ( 2007) ex- 
trapolated beyond the tabulated tem perature range of 3000 K 
to 38000 K using iBurgess & Tullvl (il992). CoUisional ioniza- 
tion rates from the ground state are from lArnaud & Rothenflugl 
(1985) and from excited states we use the general formula 
provided by Burgess & Chidichimo (1983). Autoionization is 
treated according to Arnaud & Rothenflud (Il985|) and di- 
electronic recombination according to IShull & van Steenberd 
(Il982h . The abundance is taken from Asplun d et all (l2009l) . We 
assumed partial frequency redistribution (PRD) of radiation over 
the line profiles for the H&K lines. 



3.3. Magnesium 

We use a three-level plus continuum model atom, including the 
ground state, the upper levels of the h & k lines and the contin- 
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uum. Energy levels and oscillator strengths ar e from the NIST 
database, collisional excitation rates are from Sigut & Pradhan 
([1995), collisional ionization from Seaton's semi-empirical for- 
mula (lAllenlll96 4, page 42) and photoionization cross-sections 
are from the Opacity ProjecQ. The abundance is taken from 
lAsplund et alJ (I2QQ9) . The h&k lines are computed assuming 
PRD. 



4. Cooling from lines and continue 

The basic cooling process is an electron impact excitation fol- 
lowed by a radiative deexcitation. If this photon escapes the 
atmosphere, the corresponding energy has been taken out of 
the thermal pool which corresponds to radiative cooling. In the 
corona these processes are the dominant ones: collisional excita- 
tion followed by radiative deexcitation with the photon escaping. 
Collisional de-excitation and absorption of photons (radiative 
excitation) are insignificant in comparison. This is often called 
the coronal approximation. 

In the chromosphere we have a much more complex situation 
with all rates playing a role. In strong lines most of the photons 
emitted are immediately absorbed in the same transition but the 
collisional destruction probability is often very low and the pho- 
ton may escape after a very large number of scattering events. 
We may also have photons absorbed in one line (e.g., Lyman- 
JS) escaping in another transition (e.g., H-a). A full description 
of these processes entails solving the transfer equations coupled 
with the non-LTE rate equations (or statistical equilibrium equa- 
tions if rates are assumed to be instantaneous). 

We here choose to describe the net eff'ect as a combination 
of three factors that need to be determined empirically from our 
simulations: (1) the optically thin radiative loss function (that 
gives the energy loss from the total generation of local photons 
in the absence of absorption per atom in the ionization stage 
of consideration (Hi, Can and Mgii)), (2) the probability that 
this energy escapes from the atmosphere and f the fraction of 
atoms in the given ionization stage: 



Assume a gas where 



-1 



(2) 



(3) 



Here Rji is the radiative rate coefficient from bound level j to / 
and Cjk the collisional rate coefficient from level j to k and n\ 
is the number of energy levels of the atom. These assumptions 
imply that an excited atom will nearly always fall back to its 
ground state through one or more radiative transitions. Hence, 
photon absorption will always be followed by one or more ra- 
diative deexcitations, and the gas cannot absorb energy from the 
radiation field. The only way to lose energy from the gas through 
a bound-bound transition is through a collisional excitation fol- 
lowed by one or more radiative deexcitations. All energy of the 
initial excitation is then emitted as radiation and thermal energy 
is converted into radiation. 

Therefore, as long as Eqs.[2]-[3]hold, the radiative losses Lx^ 
due to lines of atoms of species X in ionization stage m, per atom 
in ionization stage m per electron, is given by 



j-2 



Xifij 



(4) 



Qx = -LxJT)ExjT)^(T)Ax—n,p 
Nx P 



(1) 



where Qx is the radiative heating (positive Q) or cooling (nega- 
tive Q) per volume from element X (in our case H, Ca and Mg), 
Lx^(T) is the optically thin radiative loss function as function of 
temperature, T, per electron, per particle of element X in ion- 
ization stage m, Ex„^ (r) is the escape probability as function of 
some depth parameter r, ^^(T) is the fraction of element X that 

is in ionization stage m. Ax is the abundance of element X, ^ 4.1. Optically thin radiative loss function 



where A^x,„,i is the population density in the ground state of el- 
ement X, ionization stage m, xij is the energy diff'erence be- 
tween the ground state (level 1) and level j. Collisions with 
electrons are dominant under typical chromospheric conditions 
and and therefore also Lx^ will be a function of temperature 

only (apart from the factor -rp^ which is close to one). There 
are therefore reasons to believe that we may write the optically 
thin loss function as a function of temperature only. In Sec. 14.11 
we will use our numerical simulations to deduce such relations 
and compare them with the total collisional excitation from the 
ground state (Eq.|4]). 

In the chromosphere not all photons will escape the atmo- 
sphere. We parametrize this with the escape probability function 
Ex^ in Eq.[T] In Sec. 14. 21 we use our simulations to calculate this 
escape probability empirically. 

We complete the recipe by determining in Sec 14. 3 1 the factor 
^(r) in Eq.[TJ the fraction of hydrogen, calcium and magne- 
sium that resides in H i, Ca ii and Mg ii. 



is the number of hydrogen particles per gram of stellar material, 
which is a constant for a g iven set of abu ndances (=4.407x10^^ 
for the solar abundances of lAsplund et al.l (I2009D ). is the elec- 
tron number density and p is the mass density. 

The optically thin radiative loss function is equal to the col- 
lisional excitation rate in the coronal approximation and is then 
only a function of temperature. In the chromosphere this is not 
the case but the hope is that we may still write a radiative loss 
function as function of temperature only. We expect various 
routes between the levels (rather than a collisional excitation in 
one transition followed by a radiative deexcitation in the same 
transition) but as long as collisional de-excitation rates are small 
we expect the total collisional excitation from the ground state 
to be equal to the radiative losses: 



ihttp : //cdsweb . u^trasbg . f r/topbase/TheOP . html | 



From the RADYN simulation we calculate the total radiative 
losses as the sum of the net downward radiative rates multi- 
plied with the energy diff'erence of the transition, summed over 
all bound-bound transitions and the Lyman-continuum transition 
(the other continuum transitions get thin low down in the atmo- 
sphere and their influence on the energy balance can be taken 
into account in a multi-group opacity approach). 

Figure |2] shows the probability density function (PDF) of the 
total radiative losses per neutral hydrogen atom per electron as 
function of temperature for the RADYN simulatiorQ. We ex- 
clude the startup phase of the first 1000 s of the simulation and 

^ This figure, and many of the following figures, are not scatterplots 
but 2D histograms normalized for each x-axis bin. The darkness at a 
y-axis bin is thus a measure of the probability of finding an atmospheric 
point at this particular y-value given that it is in the x-axis bin. 
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Fig. 2. Probability density function of the radiative losses in 
lines and the Lyman-continuum of H i in the RADYN test at- 
mosphere as function of temperature for points above 1.7 Mm 
height. The curves give the adopted fit (red) and the total colli- 
sional excitation according to Eq.^(blue). 
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Fig. 3. Fit to the optically thin cooling per electron per hydrogen 
particle (thick solid), contributions from Lyman-a (thin black 
solid), Lyman-yS (dashed), Lyman-y (dot-dashed), Yi-a (red), 
Lyman-continuum (blue) and optically thin radiative losses from 
other elements than hydrogen (thick dotted) as function of tem- 
perature. 



only include points above 1.7 Mm height to avoid most points 
where the escape probability is small. For temperatures above 
10 kK the relation is very tight and the PDF is close to the total 
excitation curve. Below 10 kK we have a larger spread, both 
because this part includes points with a low escape probabil- 
ity and because the approximations (Eqns. [2]-[3]) start to break 
down. At low temperatures, the optical depth in the Lyman tran- 
sitions builds up quickly when integrating downwards in the at- 
mosphere and only the uppermost of these points have an escape 
probability close to one. To account for points with low escape 
probability we therefore take the upper envelope as the adopted 
optically thin radiative loss function. 

Figure [3] shows the various contributions to the radiative 
cooling in hydrogen, this time per electron per hydrogen atom 
(and not per neutral hydrogen atom) to have a normalization 
similar to the one used for coronal optically thin radiative losses. 
In most of the temperature range, the Lyman-CK transition dom- 
inates with at most 20% coming from other transitions (where 
Lyman-/? and the Lyman-continuum are the most important 
ones). At the low temperature end, the Lyman-continuum starts 
to contribute and H-a eventually dominates below 7 kK. The 
hydrogen optically thin radiation dominates over contributions 
from other elements below 32 kK. 

Figures |4]-[5] show the probability density function of the to- 
tal radiative losses per singly ionized calcium and magnesium 
atom, respectively, as function of temperature for the Bifrost 
simulation. This cooling was computed in detail using the ra- 
diative transfer code MultiSd (iLeenaarts & Carlsson! 2009) as- 
suming statistical equilibrium and plane-parallel geometry with 
each column in the snapshot treated as an independent ID atmo- 
sphere. 

The probability density function is close to the total colli- 
sional excitation curve for both Ca ii and Mg ii. In the case of 
Can this may be a bit surprising at first glance. The Can 3d^D 
doublet is metastable and we do not include transitions from it 
to the ground state. This implies Eq.[2]does not hold. However, 
for typical chromospheric electron densities, gas temperatures 



and radiation temperatures in the infrared triplet the following 
relation holds: 

Rij » ^ Cik, j > i. (5) 

This means that an atom in a Ca ii 3d ^D state is much more likely 
to get radiatively excited to one of the Ca ii 4p ^P states than to 
move to another level through collisions. Therefore, transition 
chains starting with collisional excitation from the ground state 
to a metastable level typically go through a number of radia- 
tive transitions, and ultimately end in the ground state through 
a spontaneous emission from a Can 4p^P state. The original 
collision energy is then lost from the gas and Eq. |4] still holds. 
However, net collisional rates between the excited levels give a 
larger spread and also radiative losses that sometimes are smaller 
than the total collisional excitation from the ground state. In con- 
trast to the case of hydrogen, the spread is thus not caused by 
optical depth eff'ects and hence we take the average of the values 
for our fits rather than the upper envelope. At the lowest temper- 
atures we also get points with radiative heating (not visible in 
the logarithmic plots) and the fits make a sharp dip downwards. 
See Sec.[5]for a discussion of radiative heating. 

4.2. Escape probability 

The assumption that all collisional excitations lead to emit- 
ted photons breaks down in the chromosphere, where increased 
mass density leads to increased probability of absorption of pho- 
tons and collisional deexcitation. Once radiative excitation pro- 
cesses come into play we have many more possible routes in 
the statistical equilibrium. To fully account for these processes 
we have to solve the coupled radiative transfer equations and 
non-LTE rate equations. Instead, we calculate an empirical es- 
cape probability from the ratio of the actual cooling as obtained 
from a detailed radiative transfer computation and the optically 
thin cooling function determined in Sec. 14. 1[ Ideally we would 
like to tabulate this escape probability as a function of optical 
depth in the lines. Since we are dealing with the total cooling 
over a number of lines, there is no single optical depth scale to 
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Fig. 4. Probability density function of the radiative losses in 
lines of Can as function of temperature in the Bifrost test at- 
mosphere for points above 1.3 Mm height. The curves give the 
adopted fit (red) and the total collisional excitation according to 
EqMiblue). 
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Fig. 5. Probability density function of the radiative losses in 
lines of Mg ii as function of temperature in the Bifrost test at- 
mosphere for points above 1.3 Mm height. The curves give the 
adopted fit (red) and the total collisional excitation according to 
EqM(blue). 

be chosen. For Can and Mgii we choose to tabulate the escape 
probability as a function of column mass. The rationale behind 
this choice is that the optical depth in the resonance lines is pro- 
portional to column mass as long as the ionization stage under 
consideration is the majority species. This is the case for both 
Can and Mgii, as shown in Sec. 14.31 in the region of the atmo- 
sphere where the optical depth is significant. Another advantage 
with the choice of column mass is that it is easy to calculate 
from the fundamental variables in a snapshot of an MHD simu- 
lation. For hydrogen the situation is diff'erent. The Lyman lines 
get optically thin in the lower part of the transition region where 
hydrogen gets ionized. In the ID simulation used for calculating 
the empirical coefficients, the transition region stays at a con- 
stant column mass. On a column mass scale, the hydrogen es- 
cape probability would then be a very steep function. Instead, 
we need to use a depth- variable that is proportional to the neu- 
tral hydrogen column density. We use the total column density 
of neutral hydrogen multiplied with the constant 4x10"^"^ cm^ 
which gives a depth- variable that is close to the optical depth at 
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Fig. 6. Probability density function of the empirical escape 
probability as function of approximate optical depth at Lyman-a 
line center and the adopted fit (red). Only points with a cooling 
above 5x10^ erg s"^ g"^ are included and the startup phase of 
the simulation (first 1000 s) has been excluded. 

line center of Lyman-CK. The escape probability is high to rather 
high optical depths, because of the very large scattering proba- 
bility in the Lyman transitions. 

The escape probability depends in a complicated way on the 
density, electron density, temperature and local radiation field. 
We expect it to be close to unity at low column mass density, 
and to go to zero at large column mass density. 

In order to derive an empirical relation we calculate the aver- 
age (over time for the RADYN simulation and over space for the 
Bifrost snapshot) actual cooling and the average cooling accord- 
ing to the optically thin cooling function determined in Sec. 14. II 
as function of the chosen depth-scale. We adopt this ratio as the 
empirical escape probability as function of depth. 

Figure[6] shows the probability density function of the empir- 
ical escape probability as function of approximate optical depth 
at Lyman-Qf line center together with the adopted relation. Figs. [7] 
and [8] show the probability density function of the empirical es- 
cape probability as function of column mass for Ca ii and Mg ii 
together with the adopted relation. 

Note that the adopted relations are not the averages of the 
points shown in Figs. [6]-[8] but the ratio between the average of 
the actual cooling and the average of the optically thin cool- 
ing such that points with large cooling get larger weight. The 
adopted fit may therefore look like a poor representation of the 
PDFs of Figs. [IHHl If these figures had only included the points 
with the largest cooling, the correspondence would have been 
more obvious. 

4.3. Ionization fraction 

In Sections. 14. 1144.21 we derived expressions for the radiative 
losses per atom in the relevant ionization stage per electron. In 
order to convert these losses to actual losses per volume, the frac- 
tion of atoms in the ionization stage under study (H i, Ca ii and 
Mg ii) must be known. In general this quantity is a function of 
the electron density, temperature, radiation field, and, in case of 
small transition rates, the history of the atmosphere. 

Neither Saha ionization equilibrium nor coronal equilibrium 
are valid in the chromosphere. Therefore no quick physics-based 
recipe exists that gives the ionization fraction. Fortunately, in our 
ID and 2D test atmospheres, the ionization fraction as function 
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Fig. 7. Probability density function of the empirical escape 
probability as function of column mass for Ca ii and the adopted 
fit (red). Only points with a a cooling above 5x10^ erg s"^ g"^ 
are included. 
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Fig. 8. Probability density function of the empirical escape 
probability as function of column mass for Mg ii and the adopted 
fit (red). Only points with a a cooling above 5x10^ erg s"^ g"^ 
are included. 



of temperature shows a rather small spread around an average 
value and an empirical fit to these values gives a reasonable ap- 
proximation. We also computed fits as function of temperature 
and electron density. As the ionization fraction is relatively in- 
sensitive to the electron density these two-parameter fits do not 
yield a significantly better approximation and we settle for rela- 
tions as function of temperature only. 

Figure [9] shows the probability density function of the frac- 
tion of H I, as function of temperature. There is a certain spread 
around the adopted fit which reflects the fact that the hydrogen 
ionzation is also a function of the history of the atmosphere. It is 
also clear that the coronal approximation gives a rather poor fit 
to the simulation data for temperatures below 30 kK. 

Figure [TOl shows the same relation for Ca ii. The points are 
rather close to the fitting function except for some points be- 
low 10 kK where an unusually strong radiation field may lead 
to a lower fraction of Can for some columns of the atmosphere. 
The best fit curve lies between the coronal approximation curve 
for a six-level atom (green curve) and the one for a two-level 
atom (blue curve). The reason the coronal approximation for the 
six-level atom fails is that an overpopulation of the meta- stable 




Temperature [kK] 

Fig. 9. Probability density function of the fraction of neutral H 
atoms as function of temperature in the RADYN test atmosphere 
between heights of 0.5 Mm and 2.0 Mm after the initial 1000 s. 
Curves show the adopted fit to the PDF (red), and the coronal 
approximation (blue). 
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Fig. 10. Probability density function of the fraction of Ca atoms 
in the form of Ca ii as function of temperature in the Bifrost test 
atmosphere between heights of 0.5 Mm and 2.0 Mm. Curves 
show the adopted fit to the PDF (red), the coronal approximation 
(green) and the coronal approximation with a two-level atom 
(blue). 



levels drives an upward radiative net rate in the infrared triplet 
violating the approximation of no upward radiative rates in the 
coronal approximation. In the limit where this upward rate is 
much larger than the downward rate we recover the two-level 
approximation. 

Figure [TT] shows the ionization balance of Mgii. The points 
are very close to the fitting curve over the full range of the sim- 
ulation. 



5. Radiative heating 

In the preceding section we deduced recipes to describe the ra- 
diative cooling in transitions of H i, Ca ii and Mg ii. In the cool 
phases of the chromosphere the atmosphere can be heated by 
radiation in the same transitions. 
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Fig. 11. Probability density function of the fraction of Mg atoms 
in the form of Mg ii as function of temperature in the Bifrost test 
atmosphere between heights of 0.5 Mm and 2.0 Mm. Curves 
show the adopted fit to the PDF (red), and the coronal approxi- 
mation (blue). 



Fig. 12. Probability density function of the radiative losses in 
lines of Ca ii divided by the empirically determined escape prob- 
ability as function of temperature in the Bifrost test atmosphere 
for points above 0.3 Mm height. The red curve gives the adopted 
fit. 



The total radiative heating rate (radiative flux divergence) is 
given by 



Q 



I XviJy 

Jo 



Sy)dy, 



(6) 



where Xv is the opacity, Jy is the mean intensity and Sy is 
the source function. Assuming a two-level atom and neglecting 
stimulated emission we have the following relations 



Xv = A/x„,o^o0y = — — aQ^y——Ax — p 
A^x,„ A^x P 

Sy = eBy + (\-e)Jy 



(7) ^ 



(8) 



(9) 



where A^x,„/ is the population density of the lower level /, is a 
constant, (py is the normalized absorption profile of the line, e is 
the collisional destruction probability. By is the Planck function, 
Cji is the collisional deexcitation probability, Aji is the Einstein 
coefficient for spontaneous emission, f(T) is a function of tem- 
perature and the other terms have the same meaning as in Eq. [T] 
Inserting into Eq.[6]we get 



Q = f(T)ao 



f 

Jo 



A^x A^H 
(py(Jy - By)dy—^Ax — riep. 

Nx P 



(10) 



This equation resembles Eq. [H ^-^^ is a function of temperature 
in LTE, Jy is constant with height in the optically thin regime 
such that ^ 0v(/y - By)dy is a function of temperature (but also 
a function of the source function where Jy thermalizes). At larger 
optical depths, Jy approaches By and the integral behaves simi- 
larly to the escape probability E(t) in Eq. [T] Although there are 
many approximations involved, the similarity to Eq. [T] leads us 
to adopt the same functional form for the heating as we have 
done for the cooling. We thus just use the previous fits also when 
Lx(T) is negative (corresponding to heating) and adopt for sim- 
plicity the same escape probability function. 

To test the applicability of this scheme we plotted the proba- 
bility density function of the radiative losses (positive) and gains 
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Fig. 13. Probability density function of the radiative losses in 
lines of Mgii divided by the empirically determined escape 
probability as function of temperature in the Bifrost test atmo- 
sphere for points above 0.3 Mm height. The red curve gives the 
adopted fit. 

(negative) in the Bifrost test atmosphere, divided by the empiri- 
cally defined escape probability, as function of temperature. The 
results are shown in Fig. [12] for Call and Fig. [13] for Mgll. The 
figures show that there is a fairly large spread in the losses/gains 
at low temperatures. For Can we adopt a fit that gives gains be- 
low ~ 4 kK, for Mg ii we choose zero gains/losses below ~ 5 kK. 

Heating in hydrogen transitions is mainly due to absorp- 
tion of Lyman-a photons produced by a nearby source (a strong 
shock or the transition region). This is illustrated in Fig. [141 
where there is heating due to Ly-a photons caused by the emis- 
sion at log rric = -5.2. Such behaviour cannot be modeled with 
a simple local equation as Eq. [TO] so we choose to set the losses 
and gains for hydrogen to zero below ~ 4 kK. 



6. Comparison between detailed solution and 
recipes 

We have now derived fits for the three factors that enter into 
the radiative cooling: the optically thin radiative loss function. 
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Fig. 14. Cooling as function of column mass from hydrogen 
transitions at t=3170 s in the RADYN simulation. Total cool- 
ing from the detailed calculation (black solid), from the recipes 
(black dashed), Lyman-CK (blue), H-a (red) and the Lyman con- 
tinuum (green). 



the escape probability and the ionization fraction. To test the ac- 
curacy with which we can approximate the cooling calculated 
from the detailed solution of the non-LTE rate equations we here 
compare the detailed solution with our recipe. Since the detailed 
solution from a given test atmosphere has been used to derive 
the fits we here use diff'erent snapshots for this comparison (a 
diff'erent timeseries calculated with RADYN for hydrogen and a 
snapshot at a much later time calculated with Bifrost for calcium 
and magnesium). 

Figure [T4l shows the cooling (-Q) as a function of column 
mass for hydrogen at a time in the RADYN simulation when 
there is a shock in the chromosphere. The detailed cooling is 
dominated by Lyman-a in the shock but with a significant con- 
tribution also from H-a. Behind the shock (larger column mass), 
H-Qf dominates the cooling. The total cooling is well approxi- 
mated by our recipe. In front of the shock (smaller column mass) 
there is heating by absorption of Lyman-of photons. This heating 
is not accounted for in the recipe. 

Figures [15] and [161 show the average cooling as function of 
height in the Bifrost test atmosphere. Since hydrogen cooling 
cannot be properly calculated in the Bifrost atmosphere we only 
show the recipe result for hydrogen (and use these values also 
for the calculation of the total detailed cooling) in order to show 
the relative importance of the three elements. It is clear that the 
recipes reproduce the average cooling very well. Hydrogen dom- 
inates in the upper chromosphere and transition region while 
magnesium dominates the cooling in the mid and lower chro- 
mosphere with equally strong calcium cooling at some heights. 
Note that this average cooling varies very much from snapshot 
to snapshot depending on the location of shocks in the atmo- 
sphere. These figures should thus only be used to compare the 
detailed solution with the approximate one and not as pictures 
of the mean chromospheric energy balance. 



7. Heating from incident radiation field 

Part of the optically thin radiative losses from the corona escapes 
outwards, but an equal amount of energy is emitted towards the 
sun and is absorbed in the chromosphere where it contributes to 
radiative heating. Most of this radiation is absorbed in the con- 
tinua of neutral helium and neutral hydrogen. If one assumes 



Fig. 15. Average cooling as function of height in the Bifrost test 
atmosphere. Total cooling (black), Hi (turquoise), Mgii (blue) 
and Ca ii (red). According to the detailed calculations (solid) and 
according to the recipes (dashed). 
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Fig. 16. As Fig.[T5]but for the low-mid chromosphere. 



the coronal ionization equilibrium is only dependent on temper- 
ature, then the frequency-integrated coronal losses per volume 
(positive Q means heating, as before) are given by 



2cor = -nu nJ(T). 



(11) 



Here f(T) is a positive function dependent on temperature only 
that can be pre-computed and tabulated from atomic data. The 
emissivity associated with this loss is 



7]=- 



An 



(12) 



Because the coronal radiative losses are integrated over fre- 
quency one has to choose a single representative opacity to com- 
pute the absorption. A decent choice is to use the opacity at the 
ionization edge of the ground state of neutral helium (denoted 
by a). The opacity;^ is then given by 



P A^He 



(13) 



where Nnei is the number density of neutral helium, Nue is the 
number density of helium, Aue is the abundance of helium rel- 
ative to hydrogen, Nh is the total number density of hydrogen 



particles and — is a constant dependent on the abundances. The 
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quantity Nuqi/Nuq can be pre-computed from a calibration com- 
putation as is done for H , Ca and Mg in Sec. 14.31 

The most accurate method for computing the absorption of 
coronal radiation in the chromosphere is to compute the 3D radi- 
ation field resulting from the coronal emissivity and the absorp- 
tion using the radiative transfer equation: 

dl 

^^=V-Xl (14) 

with / the intensity and s the geometrical path length along a ray. 
The solution can be computed using standard l ong or short char- 
acteristics techniques in d ecomposed domains (iHeinemann et al.l 
I2QQ6I: iHayek et al.ll201Qh . Once the intensity is known the heat- 
ing rate can be directly computed as 

Qc^=X r /dQ- 77 = 47r;t/3D-77, (15) 

with Q the solid angle. 

The method described above is accurate, but slow. A faster 
method can be obtained by treating each column in the MHD 
simulation as a plane-parallel atmosphere and assuming that at 
each location either or;^ is zero. For each column all quantities 
then become only a function of height z and Eq.[T5]reduces to 

echrfe) = Inxiz) /fcyu)dyu = AnxJiY,. (16) 

with yu the cosine of the angle with respect to the vertical. 

If we furthermore assume the emitting region is located on 
top of the absorbing region, we can ignore rays that point out- 
ward from the interior of the star in the integral of Eq. [161 The 
intensity that impinges on the top of the absorbing region for in- 
ward pointing rays is simply the integral of the emissivity along 
the z-axis from the top (zt) to the bottom (zb) of the compu- 
tational domain, weighted with the inverse of jd to account for 
slanted rays 

/(//)cor= n r n^^)^^' (17) 

Now that we have computed the incoming intensity that im- 
pinges on the top of the absorbing region we can compute the 
intensity at all heights from the solution to the transfer equation 
with zero emissivity: 

/fcyU) = /eorOu)exp-^(^>/l^'l. (18) 

The vertical optical depth r is computed from the top of the com- 
putational domain: 

T(z)= ^V(^Odz^ (19) 

The integrals in Eqs.[T7]and[T9]can both be computed efliiciently 
in parallel and need to be computed only once, independent of 
the number of yu- values in the angle discretization. 

The method gives inaccurate results if the assumption of 
coronal emission above chromospheric absorption fails. This 
happens for example when a cool finger of chromospheric gas 
protrudes slanted into the corona. 

Figure [T7] shows a comparison of the 3D and the fast ID 
method. The top panel shows the coronal emissivity 77. It is 
mostly smooth with occasional very large peaks close to the 
transition region (red curve). The emissivity for gas colder than 
30 kK (below the red curve) is negligible. 



The second panel shows the angle- averaged radiation field 
computed in 3D (/3D in Eq. (TS]). The radiation field is fairly 
smooth as it is determined by an average of the emissivity 
throughout the corona. 

The middle panel shows the angle-averaged radiation field 
treating each column as a plane-parallel atmosphere (/id in 
Eq.[T6l). Here it is far from smooth. The radiation field is high in 
columns with a large integrated emissivity, as indicated by the 
green and red vertical stripes. Note that there is cool gas located 
above coronal gas around (x, z) = (8, 3) Mm. Here the radiation 
field is highest at the top and decreases downward, while the 
emissivity panel shows that the emission is mainly coming from 
the coronal gas located below the cool finger. 

The fourth panel shows the 3D heating rate. The heating rate 
is non-zero only below the red curve, validating our assumption 
for the fast ID method of separate absorbing and emitting re- 
gions. It is fairly smooth, with quite some heating in the cool 
bubble around (x, z) = (4, 4) Mm, caused by coronal radiation 
coming in from the sides. The bottom panel shows the ID heat- 
ing rate. It shows large extrema in columns with large integrated 
emissivity, and is less smooth than the 3D heating. Because it 
is computed column-by-column, it cannot include the sideways 
radiation in the cool bubble as indicated by the low heating rate 
there compared to the 3D case. However, its overall appearance 
is quite similar to the 3D heating rate. 

Figure [18] compares the horizontally averaged heating rate 
as function of height computed from the ID recipe with the 3D 
heating rate. The recipe is always within 15% of the 3D heat- 
ing. Part of this diff'erence is caused by numerical errors because 
the strong coronal emission close to the transition region is only 
barely resolved on the grid used in the simulation. 

8. Conclusions 

We have shown that the radiative cooling in the solar chromo- 
sphere by hydrogen lines, the Lyman continuum, the H and K 
and infrared triplet lines of singly ionized calcium, and the h and 
k lines of singly ionized magnesium can be expressed fairly ac- 
curately as easily computed functions. These functions are the 
product of an optically thin emission (dependent on tempera- 
ture), an escape probability (dependent on column mass) and 
an ionization fraction (dependent on temperature). While at any 
given location in the atmosphere our recipe might deviate from 
the true cooling, it is nevertheless remarkably accurate when av- 
eraged over time and/or space as indicated by Figs.[T4l-[T6l 

We have shown that radiative heating can be approximated 
with the same functional form as radiative cooling. However, the 
heating is strongly dependent on the details of the local radiation 
field and absorption profile, so the recipe is not nearly as accu- 
rate for heating as it is for cooling. Our detailed computations 
show that radiative heating in the h and k lines of Mg 11 is negli- 
gible, and small, but present, for the lines of Ca 11. Heating due 
to hydrogen is mainly through absorption of Lyman-of photons 
next to the transition region or shock waves. This eff'ect cannot 
be modeled with our recipes. 

In addition we have presented a simple method to com- 
pute chromospheric heating by absorption of coronal radiation. 
This method approximates the diff'use component of the heat- 
ing, but yields localized spots of spurious intense heating ow- 
ing to the column-by-column approach that prevents spreading 
of radiation. It reproduces the spatially averaged heating very 
well. Given its simplicity and speed it is a good alternative to the 
slower, more accurate method that requires a formal solution of 
the transfer equation. 
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Fig. 17. Comparison of the ID and 3D method for computing 
chromospheric heating due to absorption of coronal radiation. 
All panels show a vertical slice through a snapshot of a three- 
dimensional radiation-MHD simulation performed with Bi frost. 
The red curve is located at a gas temperature of 30,000 K and in- 
dicates the location of the transition zone. All brightness ranges 
have been clipped to bring out variation at low values of the 
displayed quantities. Panels, from top to bottom: coronal emis- 
sivity 7/; angle- averaged radiation field computed in 3D; angle- 
averaged radiation field computed in ID: heating rate per mass 
with 3D radiation: heating rate per mass with ID radiation. Note 
that the 3D radiation field is computed using the full 3D volume, 
not only this 2D slice. 

The recipes have been developed for solar chromospheric 
conditions. For sunspot conditions or for other stars it would be 
necessary to redo the semi-empirical fits using detailed radiative 
transfer calculations applied to appropriate simulations. 
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Fig. 18. Horizontally-averaged radiative heating by absorption 
of coronal radiation as function of height in the 3D Bifrost test 
atmosphere from the detailed 3D computation (solid) and from 
the ID recipe (dashed). 
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